Gráfico de una Serie de tiempo

Vamos a analizar de forma descriptiva algunas serie de tiempo. Empezaremos por la serie de desempleo de los Estados Unidos que viene en el paquete TSstudio.

library(TSstudio)
data(USUnRate)
ts_info(USUnRate)
##  The USUnRate series is a ts object with 1 variable and 864 observations
##  Frequency: 12 
##  Start time: 1948 1 
##  End time: 2019 12
class(USUnRate)
## [1] "ts"
plot(USUnRate,main = "US Monthly Unemployment Rate",ylab="Unemployment Rate (%)")

La serie es medida mensual, es decir, presenta una frecuencia de 12. Qué características podemos observar? * Tendencia? * Estacionalidad? * Cíclos? * Varianza marginal no constante? Vamos a seleccionar un periodo de tiempo mas corto

unemployment <- window(USUnRate, start = c(1990,1))
ts_plot(unemployment,
          title = "US Monthly Unemployment Rate",
           Ytitle = "Unemployment Rate (%)",
           Xtitle = "Year",
          Xgrid = TRUE,
Ygrid = TRUE)

Note que aquí podemos ver varias características: Estacionalidad(no es tan evidente), tres periodos de ciclos, el primero de 1990 a 2000,el segundo de 2000 a 2007, y el tercero de 2007 a 2019. No parece tener una heterocedasticidad marginal.

Veamos ahora la tasa de desempleo de Colombia. Hay que hacerle un ajuste a la base de datos porque está en orden descendente en el tiempo.

Desempleo y Empleo Colombia

library(readxl)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
DesempleoyEmpleo <- read_excel("DesempleoyEmpleo.xlsx", range="A9:C249")
str(DesempleoyEmpleo)
## tibble [240 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Fecha          : chr [1:240] "2020-12" "2020-11" "2020-10" "2020-09" ...
##  $ Tasadeempleo   : num [1:240] 53.4 53.2 53.2 50.6 49.3 ...
##  $ Tasadedesempleo: num [1:240] 13.4 13.3 14.7 15.8 16.8 ...
DesempleoyEmpleo_1=DesempleoyEmpleo %>% map_df(rev)
tail(DesempleoyEmpleo)
## # A tibble: 6 × 3
##   Fecha   Tasadeempleo Tasadedesempleo
##   <chr>          <dbl>           <dbl>
## 1 2001-06         51.8            15.2
## 2 2001-05         51.2            14.2
## 3 2001-04         51.5            14.6
## 4 2001-03         53.0            15.7
## 5 2001-02         52.7            17.3
## 6 2001-01         53.0            16.7
head(DesempleoyEmpleo_1)
## # A tibble: 6 × 3
##   Fecha   Tasadeempleo Tasadedesempleo
##   <chr>          <dbl>           <dbl>
## 1 2001-01         53.0            16.7
## 2 2001-02         52.7            17.3
## 3 2001-03         53.0            15.7
## 4 2001-04         51.5            14.6
## 5 2001-05         51.2            14.2
## 6 2001-06         51.8            15.2
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(xts)
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
Fechas=as.yearmon(DesempleoyEmpleo_1$Fecha)
Desempleo_Col_xts=xts(x = DesempleoyEmpleo_1$Tasadedesempleo,frequency = 12,order.by = Fechas)
ts_info(Desempleo_Col_xts)
##  The Desempleo_Col_xts series is a xts object with 1 variable and 240 observations
##  Frequency: monthly 
##  Start time: Jan 2001 
##  End time: Dec 2020
plot(Desempleo_Col_xts)

ts_plot(Desempleo_Col_xts,
           title = "Tasa de Desemplo Mensual Colombia",
           Ytitle = "Tasa de Desempleo(%)",
           Xtitle = "Año",
           Xgrid = TRUE,
Ygrid = TRUE)

Qué características presenta esta serie?

Otros objetos

Note que también se puede crear un objeto tsibble y usar el paquete feast, fable y fabletools(tsibble) y timetk(tibble), lo cuales funcionan con el paquete tidyverse.

require(feasts)
## Loading required package: feasts
## Loading required package: fabletools
require(fable)
## Loading required package: fable
require(timetk)
## Loading required package: timetk
require(tsibble)
## Loading required package: tsibble
## 
## Attaching package: 'tsibble'
## The following object is masked from 'package:zoo':
## 
##     index
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
###Creación objeto tssible a partir de un objeto tibble
df_desempleo=data.frame(Desempleo=DesempleoyEmpleo_1$Tasadedesempleo,Fecha=Fechas)
tbl_desempleo=tibble(df_desempleo)
tbl_desempleo_format_fecha=tbl_desempleo
tbl_desempleo_format_fecha$Fecha=yearmonth(tbl_desempleo_format_fecha$Fecha)
###El tipo de fechas debe ser alguno que reconozca tsibble

tsbl_desempleo=as_tsibble(tbl_desempleo_format_fecha,index=Fecha)   ####La fecha en tsibble es importante

##Gráfica de tsibble
autoplot(tsbl_desempleo,Desempleo)+labs(tittle="Serie de Desempleo Colombia Mensual",y="Tasa de Deszempleo")

###Gráfica timetk
tbl_desempleo%>%plot_time_series(.value=Desempleo,.date_var=Fecha)

Análisis de Tendencias

Vamos a ver la forma de estimar la tendencia y/o eliminarla.

library(astsa)
library(TSstudio)
data(chicken)
ts_info(chicken)
##  The chicken series is a ts object with 1 variable and 180 observations
##  Frequency: 12 
##  Start time: 2001 8 
##  End time: 2016 7
plot(chicken,main="Precio Mensual de la Libra de Pollo en Estados Unidos", ylab="Precio en Centavos de Dólar")

#ts_plot(chicken)

Al parecer la serie de precios mensuales del pollo presenta una tendencia creciente al parecer lineal, es decir \[y_t=\mu_t+a_t\] o mas específicamente

\[y_t=\beta_0+\beta_1 t +a_t\]

summary(fit <- lm(chicken~time(chicken), na.action=NULL))
## 
## Call:
## lm(formula = chicken ~ time(chicken), na.action = NULL)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7411 -3.4730  0.8251  2.7738 11.5804 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.131e+03  1.624e+02  -43.91   <2e-16 ***
## time(chicken)  3.592e+00  8.084e-02   44.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.696 on 178 degrees of freedom
## Multiple R-squared:  0.9173, Adjusted R-squared:  0.9168 
## F-statistic:  1974 on 1 and 178 DF,  p-value: < 2.2e-16
plot(chicken, ylab="centavos por libra") 
abline(fit,col = "red") # Se añade la recta ajusta

###Eliminamos la tendencia con la predicción la recta
ElimiTendchick=chicken-predict(fit)
plot(ElimiTendchick,main="Serie Chicken Sin tendencia")

library(tidyverse)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:tsibble':
## 
##     interval
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(timetk)
library(tsibble)

# Se configura para gráficas plotly(# FALSE retorna ggplots y no plotly)
interactive <- FALSE
indice_chicken=as.Date(as.yearmon(tk_index(chicken)))
## Otra forma de extraer el indice estimetk::tk_index(chicken)
df_chicken=data.frame(Fecha=indice_chicken,Pollo=as.matrix(chicken))
str(df_chicken)
## 'data.frame':    180 obs. of  2 variables:
##  $ Fecha: Date, format: "2001-08-01" "2001-09-01" ...
##  $ Pollo: num  65.6 66.5 65.7 64.3 63.2 ...
tibble_chicken=tibble(df_chicken)
duplicates(tibble_chicken, key = NULL, index=Fecha)   ##Mirar si hay registros duplicados
## # A tibble: 0 × 2
## # … with 2 variables: Fecha <date>, Pollo <dbl>
## # ℹ Use `colnames()` to see all variable names
tibble_chicken_fechas_correct=tibble_chicken
tibble_chicken_fechas_correct$Fecha=yearmonth(tibble_chicken_fechas_correct$Fecha)
print(duplicates(tibble_chicken_fechas_correct, key = NULL, index=Fecha))
## # A tibble: 0 × 2
## # … with 2 variables: Fecha <mth>, Pollo <dbl>
## # ℹ Use `colnames()` to see all variable names
tsibble_chicken=tsibble(tibble_chicken_fechas_correct,index=Fecha)

tsibble_chicken
## # A tsibble: 180 x 2 [1M]
##       Fecha Pollo
##       <mth> <dbl>
##  1 2001 Aug  65.6
##  2 2001 Sep  66.5
##  3 2001 Oct  65.7
##  4 2001 Nov  64.3
##  5 2001 Dec  63.2
##  6 2002 Jan  62.9
##  7 2002 Feb  62.9
##  8 2002 Mar  62.7
##  9 2002 Apr  62.5
## 10 2002 May  63.4
## # … with 170 more rows
## # ℹ Use `print(n = ...)` to see more rows

Note que se ajusta una tendencia que se ve que no es lineal.

tibble_chicken%>%timetk::plot_time_series(Fecha, Pollo, 
                   .interactive = interactive,
                   .plotly_slider = TRUE)

###Usa Loess para hacer el ajuste de la tendencia, es decir usar smooth_vec() como versión simplificada de stats::loess()
tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.75, degree = 2))
## # A tibble: 180 × 3
##    Fecha      Pollo Pollo_ajus
##    <date>     <dbl>      <dbl>
##  1 2001-08-01  65.6       62.8
##  2 2001-09-01  66.5       63.0
##  3 2001-10-01  65.7       63.2
##  4 2001-11-01  64.3       63.4
##  5 2001-12-01  63.2       63.6
##  6 2002-01-01  62.9       63.8
##  7 2002-02-01  62.9       64.0
##  8 2002-03-01  62.7       64.3
##  9 2002-04-01  62.5       64.5
## 10 2002-05-01  63.4       64.7
## # … with 170 more rows
## # ℹ Use `print(n = ...)` to see more rows
tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.5, degree = 1))%>%
  ggplot(aes(Fecha, Pollo)) +
    geom_line() +
    geom_line(aes(y = Pollo_ajus), color = "red")

Descomposición

Vamos a usar la descomposición por medio de filtros de promedios móviles

chicken_decompo=decompose(chicken)
plot(chicken_decompo)

chicken_decompo$trend
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2001                                                                      
## 2002        NA  63.91958  63.75667  63.54250  63.32875  63.15500  63.05458
## 2003  63.64792  63.96792  64.36708  64.81542  65.31667  65.89792  66.51458
## 2004  71.92417  72.96750  73.82292  74.50458  75.07875  75.53583  75.88917
## 2005  75.35458  74.87458  74.52875  74.33750  74.18583  74.00208  73.75333
## 2006  71.02208  70.63875  70.27000  69.88542  69.53458  69.30333  69.28708
## 2007  73.70708  74.62875  75.53333  76.40667  77.19292  77.87083  78.43000
## 2008  80.89500  81.48792  82.07125  82.68125  83.38750  84.19292  85.03333
## 2009  87.24208  87.18625  86.97083  86.62875  86.23333  85.83042  85.45208
## 2010  84.68500  84.69750  84.85958  85.14083  85.44125  85.71333  85.92833
## 2011  86.39167  86.38500  86.45042  86.59625  86.84792  87.19125  87.60042
## 2012  91.02625  91.62125  92.18625  92.74917  93.34292  93.97792  94.66958
## 2013  99.52667 100.49167 101.40917 102.23292 102.95292 103.56333 104.05833
## 2014 106.44125 106.96125 107.52875 108.20167 108.95417 109.73583 110.53667
## 2015 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## 2016 113.03167        NA        NA        NA        NA        NA        NA
##            Aug       Sep       Oct       Nov       Dec
## 2001        NA        NA        NA        NA        NA
## 2002  63.03542  63.09125  63.18125  63.27625  63.41958
## 2003  67.17167  67.90875  68.76083  69.72792  70.79583
## 2004  76.14000  76.26292  76.26458  76.13750  75.82708
## 2005  73.41375  72.99042  72.48750  71.95000  71.45333
## 2006  69.53958  70.06750  70.84500  71.77125  72.74708
## 2007  78.90083  79.32750  79.69417  80.02250  80.39333
## 2008  85.76458  86.26667  86.59333  86.87833  87.12667
## 2009  85.13500  84.92125  84.84500  84.81958  84.75667
## 2010  86.08375  86.24417  86.37750  86.42792  86.42208
## 2011  88.07750  88.61125  89.17625  89.77500  90.40333
## 2012  95.41000  96.14625  96.89542  97.70167  98.58000
## 2013 104.45875 104.79708 105.15125 105.54292 105.96083
## 2014 111.32708 112.08917 112.78208 113.39792 113.91000
## 2015 114.46792 114.28542 114.03375 113.72917 113.39000
## 2016

Simule unos datos III, crear un objeto con periodicidad s=12, y extraer la descomposición, qué puede observar?

y=ts(rnorm(1000,0,1),start=c(2000,01),frequency = 4)
plot(decompose(y))

Descomposición STL

STL son las iniciales de “Seasonal and Trend decomposition using Loess”,el cual fue desarrollado por R. B. Cleveland et al. (1990).

library(feasts)
library(fable)

tsibble_chicken<-as_tsibble(chicken)
str(tsibble_chicken)
## tbl_ts [180 × 2] (S3: tbl_ts/tbl_df/tbl/data.frame)
##  $ index: mth [1:180] 2001 Aug, 2001 Sep, 2001 Oct, 2001 Nov, 2001 Dec, 2002 Jan...
##  $ value: num [1:180] 65.6 66.5 65.7 64.3 63.2 ...
##  - attr(*, "key")= tibble [1 × 1] (S3: tbl_df/tbl/data.frame)
##   ..$ .rows: list<int> [1:1] 
##   .. ..$ : int [1:180] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..@ ptype: int(0) 
##  - attr(*, "index")= chr "index"
##   ..- attr(*, "ordered")= logi TRUE
##  - attr(*, "index2")= chr "index"
##  - attr(*, "interval")= interval [1:1] 1M
##   ..@ .regular: logi TRUE
tsibble_chicken %>%
  model(
    STL(value ~ trend() +
                   season(window = "periodic"),
    robust = TRUE)) %>%
  components() %>%
  autoplot()

Note que en ambos casos se obliga a extraer un componente estacional, sin embargo puede que está componente en verdad no exista, por eso se debe verificar que en efecto hay.

set.seed(154) 
w = rnorm(200); x = cumsum(w) 
wd = w +.2; xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="Caminata Aletoria", ylab='')
lines(x, col=4); abline(h=0, col=4, lty=2); abline(a=0, b=.2, lty=2)

La diferencia ordinaria de orden 1 es , \[\nabla^1 Y_t=(1-B)^1 Y_t=Y_t-Y_{t-1}\]

dx=diff(x)
plot.ts(dx, ylim=c(-5,5), main="Serie Diferenciada", ylab='')

par(mar = c(2,2,2,2))
fit = lm(chicken~time(chicken), na.action=NULL) # Regresión sobre el tiempo
par(mfrow=c(2,1))
plot(resid(fit), type="l", main="sin tendencia") 
plot(diff(chicken), type="l", main="Primera Diferencia") 

par(mar = c(3,2,3,2))
par(mfrow=c(3,1)) # plot ACFs
acf(chicken, 48, main="ACF Pollo")
acf(resid(fit), 48, main="ACF Sin tendencia") 
acf(diff(chicken), 48, main="ACF Primera Diferencia")

Transformación de Box-Cox para Estabilizar la Varianza Marginal

En ocasiones la serie presenta varianza marginal no constante a lo largo del tiempo, lo cual hace necesario tener en cuenta tal característica. En este caso, se siguiere hacer una transformación de potencia para estabilizar la varianza. Esta familia de transformaciones se llaman transformaciones Box-Cox.

\[ f_{\lambda}(u_{t})= \begin{cases} \lambda^{-1}(u^{\lambda}_{t}-1), & \text{si $u_{t} \geq 0$, para $\lambda>0$,}\\ \ln(u_{t}), &\text{ si $u_{t}>0$, para $\lambda=0$}. \end{cases} \]

data("AirPassengers")
plot(AirPassengers)

#####Transformación Box-Cox
library(FitAR)
## Loading required package: lattice
## Loading required package: leaps
## Loading required package: ltsa
## Loading required package: bestglm
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:FitAR':
## 
##     BoxCox
## The following object is masked from 'package:astsa':
## 
##     gas
forecast::BoxCox.lambda(AirPassengers, method = "guerrero", lower = -1, upper = 3) ###Me entrega el valor de lambda 
## [1] -0.294731
##method="loglik"
FitAR::BoxCox(AirPassengers)###Me entrega una gráfica

plot(forecast::BoxCox(AirPassengers,lambda=-0.294731))###

lAirPass=log(AirPassengers)
x11()
par(mar = c(1,1,1,1))
par(mfrow=c(2,1))
plot(AirPassengers,main="Serie de Pasajeros sin Transformar")
plot(lAirPass,main="Series con Transformación BoxCox")

##Box-Cox con timetk
timetk::box_cox_vec(AirPassengers,lambda = 'auto',silent = F)
## box_cox_vec(): Using value for lambda: -0.294715585559316
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1949 2.548484 2.561375 2.588408 2.582937 2.567506 2.593720 2.615089 2.615089
## 1950 2.555038 2.577299 2.603899 2.593720 2.575381 2.616631 2.646225 2.646225
## 1951 2.610379 2.618160 2.656279 2.636912 2.648795 2.656279 2.680103 2.680103
## 1952 2.647515 2.658701 2.673640 2.659900 2.662270 2.699009 2.709884 2.720049
## 1953 2.676904 2.676904 2.715050 2.714201 2.709007 2.720866 2.737089 2.742835
## 1954 2.685298 2.668053 2.714201 2.707236 2.713347 2.737089 2.762580 2.756933
## 1955 2.720049 2.712489 2.739270 2.740706 2.741419 2.770363 2.796341 2.787869
## 1956 2.751056 2.746317 2.771524 2.769193 2.772100 2.801088 2.818144 2.814820
## 1957 2.770363 2.761963 2.792420 2.788382 2.791921 2.821786 2.837892 2.838594
## 1958 2.784223 2.772100 2.795371 2.788382 2.795857 2.826872 2.846724 2.851232
## 1959 2.794394 2.785275 2.815241 2.810978 2.820985 2.840332 2.864126 2.867216
## 1960 2.819775 2.808794 2.820583 2.836477 2.840332 2.860370 2.883509 2.879580
##           Sep      Oct      Nov      Dec
## 1949 2.595457 2.563441 2.529834 2.561375
## 1950 2.629937 2.590196 2.552878 2.602242
## 1951 2.663443 2.635540 2.611963 2.640966
## 1952 2.690331 2.671428 2.648795 2.674735
## 1953 2.715895 2.692301 2.658701 2.682201
## 1954 2.733382 2.709007 2.684272 2.709007
## 1955 2.768604 2.744238 2.715895 2.747003
## 1956 2.791921 2.765020 2.742129 2.765020
## 1957 2.814399 2.787869 2.764414 2.782096
## 1958 2.814399 2.793903 2.767420 2.782631
## 1959 2.837187 2.815659 2.795371 2.814820
## 1960 2.852177 2.836477 2.808352 2.825716

Note que ahora usamos la misma función para verificar si en verdad la varianza fue estabilizada.

FitAR::BoxCox(lAirPass)

forecast::BoxCox.lambda(lAirPass, method = "guerrero", lower = -1, upper = 2)
## [1] -0.419081
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:tidyr':
## 
##     fill
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:VGAM':
## 
##     logit
## The following object is masked from 'package:FitAR':
## 
##     Boot
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
VGAM::yeo.johnson(AirPassengers, lambda = 0)
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1949 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
## 1950 4.753590 4.844187 4.955827 4.912655 4.836282 5.010635 5.141664 5.141664
## 1951 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
## 1952 5.147494 5.198497 5.267858 5.204007 5.214936 5.389072 5.442418 5.493061
## 1953 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
## 1954 5.323010 5.241747 5.463832 5.429346 5.459586 5.579730 5.713733 5.683580
## 1955 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
## 1956 5.652489 5.627621 5.762051 5.749393 5.765191 5.926926 6.025866 6.006353
## 1957 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## 1958 5.831882 5.765191 5.894403 5.855072 5.897154 6.077642 6.198479 6.226537
## 1959 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## 1960 6.035481 5.971262 6.040255 6.135565 6.159095 6.284134 6.434547 6.408529
##           Sep      Oct      Nov      Dec
## 1949 4.919981 4.787492 4.653960 4.779123
## 1950 5.068904 4.897840 4.744932 4.948760
## 1951 5.220356 5.093750 4.990433 5.117994
## 1952 5.347108 5.257495 5.153292 5.273000
## 1953 5.472271 5.356586 5.198497 5.308268
## 1954 5.560682 5.438079 5.318120 5.438079
## 1955 5.746203 5.616771 5.472271 5.631212
## 1956 5.874931 5.726848 5.605802 5.726848
## 1957 6.003887 5.852202 5.723585 5.820083
## 1958 6.003887 5.886104 5.739793 5.823046
## 1959 6.139885 6.011267 5.894403 6.006353
## 1960 6.232448 6.135565 5.968708 6.070738
car::yjPower(AirPassengers,lambda=0)
##   [1] 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
##   [9] 4.919981 4.787492 4.653960 4.779123 4.753590 4.844187 4.955827 4.912655
##  [17] 4.836282 5.010635 5.141664 5.141664 5.068904 4.897840 4.744932 4.948760
##  [25] 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
##  [33] 5.220356 5.093750 4.990433 5.117994 5.147494 5.198497 5.267858 5.204007
##  [41] 5.214936 5.389072 5.442418 5.493061 5.347108 5.257495 5.153292 5.273000
##  [49] 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
##  [57] 5.472271 5.356586 5.198497 5.308268 5.323010 5.241747 5.463832 5.429346
##  [65] 5.459586 5.579730 5.713733 5.683580 5.560682 5.438079 5.318120 5.438079
##  [73] 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
##  [81] 5.746203 5.616771 5.472271 5.631212 5.652489 5.627621 5.762051 5.749393
##  [89] 5.765191 5.926926 6.025866 6.006353 5.874931 5.726848 5.605802 5.726848
##  [97] 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## [105] 6.003887 5.852202 5.723585 5.820083 5.831882 5.765191 5.894403 5.855072
## [113] 5.897154 6.077642 6.198479 6.226537 6.003887 5.886104 5.739793 5.823046
## [121] 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## [129] 6.139885 6.011267 5.894403 6.006353 6.035481 5.971262 6.040255 6.135565
## [137] 6.159095 6.284134 6.434547 6.408529 6.232448 6.135565 5.968708 6.070738

Ejercicio: Use la serie varve del paquete astsa para chequear si es necesario hacer transformación Box-Cox.

Vamos a correr lo mismo pero en Python para la transformación de BoxCox

import sys
print(sys.path)
## ['', '/opt/anaconda3/envs/Python38andR/bin', '/opt/anaconda3/envs/Python38andR/lib/python38.zip', '/opt/anaconda3/envs/Python38andR/lib/python3.8', '/opt/anaconda3/envs/Python38andR/lib/python3.8/lib-dynload', '/Users/sergiocalderon/.local/lib/python3.8/site-packages', '/opt/anaconda3/envs/Python38andR/lib/python3.8/site-packages', '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/reticulate/python']
import pandas as pd
import matplotlib.pylab as plt
data = pd.read_csv('AirPassengers.csv')
print(data)
##        Month  NPassengers
## 0    1949-01          112
## 1    1949-02          118
## 2    1949-03          132
## 3    1949-04          129
## 4    1949-05          121
## ..       ...          ...
## 139  1960-08          606
## 140  1960-09          508
## 141  1960-10          461
## 142  1960-11          390
## 143  1960-12          432
## 
## [144 rows x 2 columns]
print('\n Data Types:')
## 
##  Data Types:
print(data.dtypes)
## Month          object
## NPassengers     int64
## dtype: object
con=data['Month']
data['Month']=pd.to_datetime(data['Month'])
data
##          Month  NPassengers
## 0   1949-01-01          112
## 1   1949-02-01          118
## 2   1949-03-01          132
## 3   1949-04-01          129
## 4   1949-05-01          121
## ..         ...          ...
## 139 1960-08-01          606
## 140 1960-09-01          508
## 141 1960-10-01          461
## 142 1960-11-01          390
## 143 1960-12-01          432
## 
## [144 rows x 2 columns]
pasajeros=data.set_index('Month')
#check datatype of index

#convert to time series:
ts = pasajeros['NPassengers']
ts.head(10)

####Graficar la Serie#####
#plt.plot(ts)
#plt.title('AirPassengers')
#plt.show()
## Month
## 1949-01-01    112
## 1949-02-01    118
## 1949-03-01    132
## 1949-04-01    129
## 1949-05-01    121
## 1949-06-01    135
## 1949-07-01    148
## 1949-08-01    148
## 1949-09-01    136
## 1949-10-01    119
## Name: NPassengers, dtype: int64
print(r.AirPassengers)
## [112.0, 118.0, 132.0, 129.0, 121.0, 135.0, 148.0, 148.0, 136.0, 119.0, 104.0, 118.0, 115.0, 126.0, 141.0, 135.0, 125.0, 149.0, 170.0, 170.0, 158.0, 133.0, 114.0, 140.0, 145.0, 150.0, 178.0, 163.0, 172.0, 178.0, 199.0, 199.0, 184.0, 162.0, 146.0, 166.0, 171.0, 180.0, 193.0, 181.0, 183.0, 218.0, 230.0, 242.0, 209.0, 191.0, 172.0, 194.0, 196.0, 196.0, 236.0, 235.0, 229.0, 243.0, 264.0, 272.0, 237.0, 211.0, 180.0, 201.0, 204.0, 188.0, 235.0, 227.0, 234.0, 264.0, 302.0, 293.0, 259.0, 229.0, 203.0, 229.0, 242.0, 233.0, 267.0, 269.0, 270.0, 315.0, 364.0, 347.0, 312.0, 274.0, 237.0, 278.0, 284.0, 277.0, 317.0, 313.0, 318.0, 374.0, 413.0, 405.0, 355.0, 306.0, 271.0, 306.0, 315.0, 301.0, 356.0, 348.0, 355.0, 422.0, 465.0, 467.0, 404.0, 347.0, 305.0, 336.0, 340.0, 318.0, 362.0, 348.0, 363.0, 435.0, 491.0, 505.0, 404.0, 359.0, 310.0, 337.0, 360.0, 342.0, 406.0, 396.0, 420.0, 472.0, 548.0, 559.0, 463.0, 407.0, 362.0, 405.0, 417.0, 391.0, 419.0, 461.0, 472.0, 535.0, 622.0, 606.0, 508.0, 461.0, 390.0, 432.0]
print(r.lAirPass)
## [4.718498871295094, 4.770684624465665, 4.882801922586371, 4.859812404361672, 4.795790545596741, 4.90527477843843, 4.997212273764115, 4.997212273764115, 4.912654885736052, 4.77912349311153, 4.6443908991413725, 4.770684624465665, 4.74493212836325, 4.836281906951478, 4.948759890378168, 4.90527477843843, 4.8283137373023015, 5.003946305945459, 5.135798437050262, 5.135798437050262, 5.062595033026967, 4.890349128221754, 4.736198448394496, 4.941642422609304, 4.976733742420574, 5.0106352940962555, 5.181783550292085, 5.093750200806762, 5.147494476813453, 5.181783550292085, 5.293304824724492, 5.293304824724492, 5.214935757608986, 5.087596335232384, 4.983606621708336, 5.111987788356544, 5.14166355650266, 5.19295685089021, 5.262690188904886, 5.198497031265826, 5.209486152841421, 5.384495062789089, 5.438079308923196, 5.488937726156687, 5.342334251964811, 5.25227342804663, 5.147494476813453, 5.267858159063328, 5.278114659230517, 5.278114659230517, 5.4638318050256105, 5.459585514144159, 5.43372200355424, 5.493061443340548, 5.575949103146316, 5.605802066295998, 5.4680601411351315, 5.351858133476067, 5.19295685089021, 5.303304908059076, 5.318119993844216, 5.236441962829949, 5.459585514144159, 5.424950017481403, 5.455321115357702, 5.575949103146316, 5.71042701737487, 5.680172609017068, 5.556828061699537, 5.43372200355424, 5.313205979041787, 5.43372200355424, 5.488937726156687, 5.4510384535657, 5.58724865840025, 5.594711379601839, 5.598421958998375, 5.752572638825633, 5.8971538676367405, 5.849324779946859, 5.7430031878094825, 5.6131281063880705, 5.4680601411351315, 5.627621113690637, 5.648974238161206, 5.6240175061873385, 5.75890177387728, 5.746203190540153, 5.762051382780177, 5.924255797414532, 6.023447592961033, 6.003887067106539, 5.872117789475416, 5.723585101952381, 5.602118820879701, 5.723585101952381, 5.752572638825633, 5.707110264748875, 5.87493073085203, 5.8522024797744745, 5.872117789475416, 6.045005314036012, 6.142037405587356, 6.1463292576688975, 6.0014148779611505, 5.849324779946859, 5.720311776607412, 5.817111159963204, 5.8289456176102075, 5.762051382780177, 5.8916442118257715, 5.8522024797744745, 5.8944028342648505, 6.075346031088684, 6.19644412779452, 6.22455842927536, 6.0014148779611505, 5.883322388488279, 5.736572297479192, 5.820082930352362, 5.886104031450156, 5.834810737062605, 6.0063531596017325, 5.981414211254481, 6.040254711277414, 6.156978985585555, 6.306275286948016, 6.326149473155099, 6.137727054086234, 6.008813185442595, 5.8916442118257715, 6.003887067106539, 6.0330862217988015, 5.968707559985366, 6.037870919922137, 6.133398042996649, 6.156978985585555, 6.282266746896006, 6.432940092739179, 6.406879986069314, 6.230481447578482, 6.133398042996649, 5.966146739123692, 6.068425588244111]
py$data
##                   Month NPassengers
## 1   1948-12-31 19:00:00         112
## 2   1949-01-31 19:00:00         118
## 3   1949-02-28 19:00:00         132
## 4   1949-03-31 19:00:00         129
## 5   1949-04-30 19:00:00         121
## 6   1949-05-31 19:00:00         135
## 7   1949-06-30 19:00:00         148
## 8   1949-07-31 19:00:00         148
## 9   1949-08-31 19:00:00         136
## 10  1949-09-30 19:00:00         119
## 11  1949-10-31 19:00:00         104
## 12  1949-11-30 19:00:00         118
## 13  1949-12-31 19:00:00         115
## 14  1950-01-31 19:00:00         126
## 15  1950-02-28 19:00:00         141
## 16  1950-03-31 19:00:00         135
## 17  1950-04-30 19:00:00         125
## 18  1950-05-31 19:00:00         149
## 19  1950-06-30 19:00:00         170
## 20  1950-07-31 19:00:00         170
## 21  1950-08-31 19:00:00         158
## 22  1950-09-30 19:00:00         133
## 23  1950-10-31 19:00:00         114
## 24  1950-11-30 19:00:00         140
## 25  1950-12-31 19:00:00         145
## 26  1951-01-31 19:00:00         150
## 27  1951-02-28 19:00:00         178
## 28  1951-03-31 19:00:00         163
## 29  1951-04-30 19:00:00         172
## 30  1951-05-31 19:00:00         178
## 31  1951-06-30 19:00:00         199
## 32  1951-07-31 19:00:00         199
## 33  1951-08-31 19:00:00         184
## 34  1951-09-30 19:00:00         162
## 35  1951-10-31 19:00:00         146
## 36  1951-11-30 19:00:00         166
## 37  1951-12-31 19:00:00         171
## 38  1952-01-31 19:00:00         180
## 39  1952-02-29 19:00:00         193
## 40  1952-03-31 19:00:00         181
## 41  1952-04-30 19:00:00         183
## 42  1952-05-31 19:00:00         218
## 43  1952-06-30 19:00:00         230
## 44  1952-07-31 19:00:00         242
## 45  1952-08-31 19:00:00         209
## 46  1952-09-30 19:00:00         191
## 47  1952-10-31 19:00:00         172
## 48  1952-11-30 19:00:00         194
## 49  1952-12-31 19:00:00         196
## 50  1953-01-31 19:00:00         196
## 51  1953-02-28 19:00:00         236
## 52  1953-03-31 19:00:00         235
## 53  1953-04-30 19:00:00         229
## 54  1953-05-31 19:00:00         243
## 55  1953-06-30 19:00:00         264
## 56  1953-07-31 19:00:00         272
## 57  1953-08-31 19:00:00         237
## 58  1953-09-30 19:00:00         211
## 59  1953-10-31 19:00:00         180
## 60  1953-11-30 19:00:00         201
## 61  1953-12-31 19:00:00         204
## 62  1954-01-31 19:00:00         188
## 63  1954-02-28 19:00:00         235
## 64  1954-03-31 19:00:00         227
## 65  1954-04-30 19:00:00         234
## 66  1954-05-31 19:00:00         264
## 67  1954-06-30 19:00:00         302
## 68  1954-07-31 19:00:00         293
## 69  1954-08-31 19:00:00         259
## 70  1954-09-30 19:00:00         229
## 71  1954-10-31 19:00:00         203
## 72  1954-11-30 19:00:00         229
## 73  1954-12-31 19:00:00         242
## 74  1955-01-31 19:00:00         233
## 75  1955-02-28 19:00:00         267
## 76  1955-03-31 19:00:00         269
## 77  1955-04-30 19:00:00         270
## 78  1955-05-31 19:00:00         315
## 79  1955-06-30 19:00:00         364
## 80  1955-07-31 19:00:00         347
## 81  1955-08-31 19:00:00         312
## 82  1955-09-30 19:00:00         274
## 83  1955-10-31 19:00:00         237
## 84  1955-11-30 19:00:00         278
## 85  1955-12-31 19:00:00         284
## 86  1956-01-31 19:00:00         277
## 87  1956-02-29 19:00:00         317
## 88  1956-03-31 19:00:00         313
## 89  1956-04-30 19:00:00         318
## 90  1956-05-31 19:00:00         374
## 91  1956-06-30 19:00:00         413
## 92  1956-07-31 19:00:00         405
## 93  1956-08-31 19:00:00         355
## 94  1956-09-30 19:00:00         306
## 95  1956-10-31 19:00:00         271
## 96  1956-11-30 19:00:00         306
## 97  1956-12-31 19:00:00         315
## 98  1957-01-31 19:00:00         301
## 99  1957-02-28 19:00:00         356
## 100 1957-03-31 19:00:00         348
## 101 1957-04-30 19:00:00         355
## 102 1957-05-31 19:00:00         422
## 103 1957-06-30 19:00:00         465
## 104 1957-07-31 19:00:00         467
## 105 1957-08-31 19:00:00         404
## 106 1957-09-30 19:00:00         347
## 107 1957-10-31 19:00:00         305
## 108 1957-11-30 19:00:00         336
## 109 1957-12-31 19:00:00         340
## 110 1958-01-31 19:00:00         318
## 111 1958-02-28 19:00:00         362
## 112 1958-03-31 19:00:00         348
## 113 1958-04-30 19:00:00         363
## 114 1958-05-31 19:00:00         435
## 115 1958-06-30 19:00:00         491
## 116 1958-07-31 19:00:00         505
## 117 1958-08-31 19:00:00         404
## 118 1958-09-30 19:00:00         359
## 119 1958-10-31 19:00:00         310
## 120 1958-11-30 19:00:00         337
## 121 1958-12-31 19:00:00         360
## 122 1959-01-31 19:00:00         342
## 123 1959-02-28 19:00:00         406
## 124 1959-03-31 19:00:00         396
## 125 1959-04-30 19:00:00         420
## 126 1959-05-31 19:00:00         472
## 127 1959-06-30 19:00:00         548
## 128 1959-07-31 19:00:00         559
## 129 1959-08-31 19:00:00         463
## 130 1959-09-30 19:00:00         407
## 131 1959-10-31 19:00:00         362
## 132 1959-11-30 19:00:00         405
## 133 1959-12-31 19:00:00         417
## 134 1960-01-31 19:00:00         391
## 135 1960-02-29 19:00:00         419
## 136 1960-03-31 19:00:00         461
## 137 1960-04-30 19:00:00         472
## 138 1960-05-31 19:00:00         535
## 139 1960-06-30 19:00:00         622
## 140 1960-07-31 19:00:00         606
## 141 1960-08-31 19:00:00         508
## 142 1960-09-30 19:00:00         461
## 143 1960-10-31 19:00:00         390
## 144 1960-11-30 19:00:00         432
library(ggplot2)
ggplot2::ggplot(data = py$data,aes(x=Month,y=NPassengers) )+geom_line()

Gráficas de Retardos

Vamos a hacer gráficos de dispersión para chequear que tipos de relaciones hay entre los retardos de la variable interés. Vamos a trabajar con algunas series, por ejemplo: * Indice ambiental mensual (soi)(Southern Oscillation Index), el cual mide los cambios en la presión del aire, relacionados con las temperaturas de la superficie del mar en el Océano Pacífico central. * la serie rec (reclutamiento asociada al soi), número de nuevos peces. * Consumo mensual de gas natural en EE. UU.(USgas) medido en Billones de pies cúbicos. Esto permite chequear si hay posibles relaciones no-lineales.

library(astsa)
data("soi")
ts_info(soi)
##  The soi series is a ts object with 1 variable and 453 observations
##  Frequency: 12 
##  Start time: 1950 1 
##  End time: 1987 9
?soi

par(mar = c(2,2,2,2))
plot(soi, main="Indice soi")

par(mar = c(3,2,3,2))
astsa::lag1.plot(soi, 12)  ###El 12 indica cuantos retardos y_t-k contra y_t 

###Hacer la gráfica con x11()

En el gráfico de dispersión podemos ver que se muestra un ajuste no paramétrico, de la posible relación entre las variables al igual que una estimación de la autocorrelación entre \(s_t\) y \(s_{t-h}\). Vemos que varias de las relaciones exploradas parecen ser lineales. Uno puede observar que existen relaciones lineales positivas en los rezagos \(h = 1, 2, 11, 12\), mientras que negativas en los rezagos \(h=6,7\), las demás parecen ser no significativas o no lineales.

#pdf('/Users/sergiocalderon/Documents/Documentos - iMac de Sergio/Documentos iMac Sergio/Notas de Clase/Notas de clase/Notas de Clase Series de Tiempo Univariadas/Graficas/DispersionSoiRec.pdf',paper="USr")
par(mar = c(3,2,3,2))
lag2.plot(soi, rec, 8)   #El 2 de lag2.plot es porque intervenienen dos serie de tiempo.

#dev.off()
lag2.plot(soi, rec, 8,corr=F)

Note también que el gráfico de dispersión del índice de nuevos peces con retardos del indice soi nos muestra también relaciones posiblemente lineales, al igual que no lineales. Ahora, podemos crear el gráfico de autocorrelación simple que nos permite estimar y graficar las autocorrelaciones para diferentes rezagos.

Note que con la función ts_lags de TSstudio podemos hacer una gráfica similar.

ts_lags(soi,lags=1:12)

La función de autocorrelación simple o gráfico ACF y PACF

Cuando el proceso es estacionario, o al menos no presenta tendencia, podemos usar el gráfico acf para explorar las posibles relaciones lineales a diferentes rezagos. En seguida mostramos la función de autocorrelación para el índice soi y la serie de nuevos peces.

par(mfrow=c(2,1))
par(mar = c(2.7,2,2.7,2))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")

#dev.off()
ccf(soi, rec, 48, main="SOI vs Recruitment", ylab="CCF")

#Corr(rec_{t},soi_{t-h})
#Corr(soi_{t},rec_{t-h},)
index_soi_rec=as.Date(as.yearmon(tk_index(soi)))

df_soi_rec=data.frame(Fecha_soi_rec=index_soi_rec,soi=as.matrix(soi),rec=as.matrix(rec))

tibble_soi_rec=tibble(df_soi_rec)
tsibble::duplicates(tibble_soi_rec, key = NULL, index=Fecha_soi_rec)   ##Mirar si hay registros duplicados
## # A tibble: 0 × 3
## # … with 3 variables: Fecha_soi_rec <date>, soi <dbl>, rec <dbl>
## # ℹ Use `colnames()` to see all variable names
tibble_soi_rec%>%plot_acf_diagnostics(Fecha_soi_rec,soi,.ccf_vars = rec,.lags = 36)
## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

En ambas gráficas podemos ver que se muestran periodicidades en las correlaciones que corresponden a valores separados por 12 unidades. También podemos ver que las observaciones con 12 meses o un año de diferencia están fuertemente correlacionadas positivamente, al igual que las observaciones en múltiplos como 24, 36, 48,\(\cdots\) Las observaciones separadas por seis meses están correlacionadas negativamente, lo que muestra que las excursiones positivas tienden a asociarse con las excursiones negativas a los seis meses ver Shumway2017 capítulo 1.

###AMI Del los libros H. Kantz and T. Schreiber: Nonlinear Time series Analysis (Cambridge university press) H. Abarbanel: Analysis of observed chaotic data (Springer, 1996) y NONLINEAR TIME SERIES ANALYSIS HOLGER KANTZ AND THOMAS SCHREIBER. Cambrige University Press 2003.

Ahora utilizaremos los paquetes nonlinearTseries y tseriesChaos para computar el average mutual information(AMI) o La información mutua promedio (AMI, la cual mide cuánto nos dice una variable aleatoria sobre otra, el cual se define como:

\[I(X;Y)=\sum_{i}\sum_{j}p(x_i,y_j)\log_2(\frac{p(x_i,y_j)}{p(x_i)p(y_j)}).\] En el contexto del análisis de series de tiempo, AMI ayuda a cuantificar la cantidad de conocimiento obtenido sobre el valor de \(X_{t+d}\) al observar \(X_t\). Equivalentemente, el AMI es una medida de qué tanto el conocimiento de \(X\) reduce la incertidumbre acerca de \(Y\). Esto implica que \(I(X,Y)=0\) si y sólo si \(X\) y \(Y\) son variables aletorias independientes. I(X; Y ) describe la información que la medición \(X_t\) en el tiempo \(t\) aporta a la medición \(X_{t+d}\) en el tiempo \(t + d\). Si se elige d como el valor alrededor del primer mínimo del AMI, entonces \(Y{t}\) e \(y_{t+d}\) son parcialmente pero no totalmente independientes.

Vamos a simular una serie de la forma \[x_t=\frac{x_t-1}{x_{t-3}^2+1}+\epsilon_t\] y a trabajar con la serie de linces Candienses del paquete astsa lynx.

library(nonlinearTseries)
## 
## Attaching package: 'nonlinearTseries'
## The following object is masked from 'package:fabletools':
## 
##     estimate
## The following object is masked from 'package:grDevices':
## 
##     contourLines
library(tseriesChaos)
et=rnorm(1100,0,1)
xt=rep(0,1100)
for(t in 13:1100)
  {
  xt[t]=(xt[t-1]-1)/(xt[t-12]^2+1)+et[t]
  }
xtsimul=as.ts(xt[101:1100])
length(xtsimul)
## [1] 1000
plot(xtsimul)

acf(xtsimul)

par(mar = c(3,2,3,2))
astsa::lag1.plot(xtsimul, 6) 

nonlinearTseries::mutualInformation(xtsimul,lag.max = 100,n.partitions = 50,units = "Bits",do.plot = TRUE) #c("Nats", "Bits", "Bans")

## $time.lag
##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100
## 
## $mutual.information
##   [1] 4.7211355 0.9086685 0.7868768 0.7146046 0.7234251 0.7393485 0.7257831
##   [8] 0.7204692 0.7530212 0.7388565 0.6987539 0.7426234 0.8363724 0.7918314
##  [15] 0.7495687 0.7384157 0.7227720 0.7306765 0.7185728 0.7490403 0.7888126
##  [22] 0.7335948 0.7676440 0.7241823 0.7728012 0.7690717 0.7278799 0.7436484
##  [29] 0.7654490 0.7561503 0.7324961 0.7820421 0.7413698 0.7871444 0.7796870
##  [36] 0.7356812 0.7411850 0.7501907 0.7512091 0.7691925 0.7466313 0.7425748
##  [43] 0.7152309 0.7567790 0.7658219 0.7516127 0.7371519 0.7288299 0.7622234
##  [50] 0.7844343 0.7584129 0.7388515 0.7564102 0.7690676 0.7334833 0.7605943
##  [57] 0.7808753 0.7568963 0.7571249 0.7135538 0.7925615 0.7584517 0.7768902
##  [64] 0.7298138 0.7462323 0.7553250 0.7728154 0.7850680 0.7464423 0.7543733
##  [71] 0.7465533 0.7813147 0.7620542 0.7528362 0.7663679 0.8111562 0.7237969
##  [78] 0.7386509 0.7523882 0.7398758 0.7844187 0.7326143 0.7380267 0.7455960
##  [85] 0.7691433 0.7309565 0.7479580 0.7314371 0.7583339 0.8090117 0.7704502
##  [92] 0.7519965 0.7847387 0.7292170 0.7291928 0.7169975 0.7426216 0.7759557
##  [99] 0.7898280 0.7688683 0.7700155
## 
## $units
## [1] "Bits"
## 
## $n.partitions
## [1] 50
## 
## attr(,"class")
## [1] "mutualInf"
pacf(xtsimul)

tseriesChaos::mutual(xtsimul, partitions = 50, lag.max = 100, plot=TRUE)

plot(astsa::Lynx)

acf(astsa::Lynx)

par(mar = c(3,2,3,2))
astsa::lag1.plot(astsa::Lynx, 12) 

tseriesChaos::mutual(astsa::Lynx, partitions = 50, lag.max = 10, plot=TRUE)

Detección de cíclos y estacionalidades

Podemos usar la función TSstudio::ts_heatmap para crear un mapa de calor. Este es un gráfico tridimensional, en donde en el eje y están los meses, y en el eje x están los años. Note en este caso que los meses de Diciembre,Enero, Febrero y Marzo presentan los valores mas oscuros, es decir los valores mas grandes a lo largo de los años, en contraste con los meses de Mayo a Septiembre que presentan colores mas claros. Este es un típico comportamiento de la presencia de un cíclo estacional en la serie. El flujo de color es horizontal.

TSstudio::USgas
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2000 2510.5 2330.7 2050.6 1783.3 1632.9 1513.1 1525.6 1653.1 1475.0 1567.8
## 2001 2677.0 2309.5 2246.6 1807.2 1522.4 1444.4 1598.1 1669.2 1494.1 1649.1
## 2002 2487.6 2242.4 2258.4 1881.0 1611.5 1591.4 1748.4 1725.7 1542.2 1645.9
## 2003 2700.5 2500.3 2197.9 1743.5 1514.7 1368.4 1600.5 1651.6 1428.6 1553.2
## 2004 2675.8 2511.1 2100.9 1745.2 1573.0 1483.7 1584.9 1578.0 1482.2 1557.2
## 2005 2561.9 2243.0 2205.8 1724.9 1522.6 1534.1 1686.6 1695.1 1422.5 1428.2
## 2006 2165.3 2144.4 2126.4 1681.0 1526.3 1550.9 1758.7 1751.7 1462.1 1644.2
## 2007 2475.6 2567.0 2128.8 1810.1 1559.1 1555.2 1659.9 1896.1 1590.5 1627.8
## 2008 2734.0 2503.4 2278.2 1823.9 1576.4 1604.2 1708.6 1682.9 1460.9 1635.8
## 2009 2729.7 2332.5 2170.7 1741.3 1504.0 1527.8 1658.0 1736.5 1575.0 1666.5
## 2010 2809.8 2481.0 2142.9 1691.8 1617.3 1649.5 1825.8 1878.9 1637.5 1664.9
## 2011 2888.6 2452.4 2230.5 1825.0 1667.4 1657.3 1890.5 1891.8 1655.6 1744.5
## 2012 2756.2 2500.7 2127.8 1953.1 1873.8 1868.4 2069.8 2008.8 1807.2 1901.1
## 2013 2878.8 2567.2 2521.1 1967.5 1752.5 1742.9 1926.3 1927.4 1767.0 1866.8
## 2014 3204.1 2741.2 2557.9 1961.7 1810.2 1745.4 1881.0 1933.1 1809.3 1912.8
## 2015 3115.0 2925.2 2591.3 2007.9 1858.1 1899.9 2067.7 2052.7 1901.3 1987.3
## 2016 3091.7 2652.3 2356.3 2083.8 1965.8 2000.7 2186.6 2208.4 1947.8 1925.2
## 2017 2914.2 2340.6 2523.7 1932.5 1892.5 1910.9 2142.1 2094.3 1920.9 2032.0
## 2018 3335.0 2705.9 2792.6 2346.3 2050.9 2058.7 2344.6 2307.7 2151.5 2279.1
## 2019 3399.9 2999.2 2899.9 2201.1 2121.0 2115.2 2407.5 2437.2 2215.6 2472.3
##         Nov    Dec
## 2000 1908.5 2587.5
## 2001 1701.0 2120.2
## 2002 1913.6 2378.9
## 2003 1753.6 2263.7
## 2004 1782.8 2327.7
## 2005 1663.4 2326.4
## 2006 1765.4 2122.8
## 2007 1834.5 2399.2
## 2008 1868.9 2399.7
## 2009 1776.2 2491.9
## 2010 1973.3 2714.1
## 2011 2031.9 2541.9
## 2012 2167.8 2503.9
## 2013 2316.9 2920.8
## 2014 2357.5 2679.2
## 2015 2249.1 2588.2
## 2016 2159.4 2866.3
## 2017 2357.7 3084.5
## 2018 2709.9 2993.1
## 2019
 ts_plot(USgas,
           title = "US Monthly Natural Gas consumption",
           Ytitle = "Billion Cubic Feet",
           Xtitle = "Year",
           Xgrid = TRUE,
           Ygrid = TRUE)
TSstudio::ts_heatmap(USgas,title = "Mapa de Calor - Consumo de Gas Natural en EEUU")

Voy enfocarme en un periodo de la serie de desempleo.

 unemployment <- window(USUnRate, start = c(1990,1))
   ts_plot(unemployment,
           title = "US Monthly Unemployment Rate",
           Ytitle = "Unemployment Rate (%)",
           Xtitle = "Year",
           Xgrid = TRUE,
Ygrid = TRUE)
ts_info(unemployment)
##  The unemployment series is a ts object with 1 variable and 360 observations
##  Frequency: 12 
##  Start time: 1990 1 
##  End time: 2019 12
ts_heatmap(unemployment,title = "Mapa de Calor - Tasa de Desempleo EEUU Subserie")
#ts_heatmap(USUnRate,title = "Mapa de Calor - Tasa de Desempleo EEUU")

En este ejemplo, el flujo de color de la Tasa de Desempleo es vertical, lo que indica el estado del ciclo. En este caso, las franjas verticales más claras representan el final de un ciclo y el comienzo del siguiente. Asimismo, las franjas verticales más oscuras representan los picos del ciclo.

Explorando mas herramientas para detección de Estacionalidad

Medidas descriptivas

USgas_df <- data.frame(year = floor(time(USgas)), month = cycle(USgas),USgas = as.numeric(USgas))

USgas_df$month <- factor(month.abb[USgas_df$month], levels = month.abb)

library(dplyr)
USgas_summary <- USgas_df %>%group_by(month) %>%summarise(mean= mean(USgas),sd = sd(USgas))
USgas_summary
## # A tibble: 12 × 3
##    month  mean    sd
##    <fct> <dbl> <dbl>
##  1 Jan   2806.  309.
##  2 Feb   2502.  223.
##  3 Mar   2325.  241.
##  4 Apr   1886.  175.
##  5 May   1708.  194.
##  6 Jun   1691.  216.
##  7 Jul   1864.  261.
##  8 Aug   1889.  237.
##  9 Sep   1687.  242.
## 10 Oct   1788.  260.
## 11 Nov   2015.  284.
## 12 Dec   2543.  278.
 library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
  plot_ly (data = USgas_summary, x = ~ month, y = ~ mean, type = "bar", name   = "Mean") %>%
   layout (title = "USgas - Monthly Average", yaxis =list(title = "Mean",   range = c(1500, 2700)))
  monthplot(USgas)

  ####También lo podemos hacer usando el paquete feasts y el objeto tsibble
  USgas_tsbl=as_tsibble(USgas)
  require(feasts)
  USgas_tsbl=as_tsibble(USgas)
  USgas_tsbl%>%gg_subseries(value)###Puede usar el argumento period=12 y da el mismo resultado, lo que significa es que se pueden agrupar las observaciones que están cada 12.

Note que basados en las estadísticas descriptivas podemos ver que las medias son distintas para algunos meses, incluso sin quitar la tendencia. A una misma conclusión llegamos basados en los gráficos. Estas son típicas características de que hay presente un ciclo estacional en la serie.

Explorando múltiples estacionalidades

Pueden haber múltiples estacionalidades(por lo general ocurren en series de alta frecuencia: diaria, cada hora, cada media y asi sucesivamente.) o una única estacionalidad, o inclusive cíclos que no se ven con facilidad.

Vamos a trabajar la base datos relacionada con el sistema nacional de transmisión de electricidad del Reino Unido. Más específicamente con la demanda de energía de forma horaria.

library(UKgrid)
require(TSstudio)
require(timetk)
require(feasts)
require(tsibble)
require(plotly)
UKgrid_xts <- extract_grid(type = "xts",
                              columns = "ND",
                              aggregate = "hourly",
                              na.rm = TRUE)
#extract_grid solo funciona para el conjunto de datos UKgrid
ts_plot(UKgrid_xts,
            title = "National Hourly Demand UK Grid",
            Ytitle = "Megawatts",
            Xtitle = "Year",
            Xgrid = TRUE,
            Ygrid = TRUE)
UKgrid_tstbl <- extract_grid(type = "tsibble",
                              columns = "ND",
                              aggregate = "hourly",
                              na.rm = TRUE)

UKgrid_tbl <-as_tibble(UKgrid_tstbl)

Como indicamos anteriormente, la primera indicación de la posible existencia de múltiples patrones estacionales en la serie es una frecuencia alta, como por diaria, horaria y en minutos. En esos casos, hay más de una forma de establecer la frecuencia de la serie. Por ejemplo, si capturamos una serie de tiempo de frecuencia diaria, la frecuencia de la serie se puede configurar de la siguiente manera: * Diariamente (o 365), asumiendo que el ciclo más apropiado es un año completo. * Días de semana (o 7) siempre que la oscilación del día de la semana sea más dominante que la del ciclo de año completo.

Al utilizar estadísticas descriptivas con este tipo de series, tendrá sentido aplicar este método para cada frecuencia potencial de la serie (o al menos las principales) con el fin de examinar si existe una indicación del patrón estacional.

Por ejemplo, UKgrid es una serie de tiempo por horas, que la marca automáticamente como sospechosa de tener múltiples patrones estacionales. Potencialmente, como se mencionó anteriormente, la demanda horaria de electricidad podría tener tres patrones estacionales diferentes:

  • Horaria: Este es probablemente el principal patrón estacional de la serie, ya que existe una relación directa entre la demanda de energía eléctrica y la hora del día (hay alta demanda durante el día y baja demanda durante la noche o al contrario).
  • Día de la semana: La demanda de electricidad a lo largo del día se deriva, potencialmente, del día de la semana. Tendría sentido esperar un alto consumo durante los días laborables y una menor tasa de consumo durante el fin de semana.
  • Mensual: como los patrones climáticos varían a lo largo del año, la cantidad de luz del día y otros factores estacionales podrían afectar la demanda de electricidad.

Usando el paquete lubridate crearemos las características.

library(xts)
UKgrid_df <- data.frame(time = zoo::index(UKgrid_xts), UKgrid=as.numeric(UKgrid_xts))
str(UKgrid_df)
## 'data.frame':    127296 obs. of  2 variables:
##  $ time  : POSIXct, format: "2005-04-01 00:00:00" "2005-04-01 01:00:00" ...
##  $ UKgrid: num  65080 68207 69172 66769 64660 ...

Ahora crearemos características estacionales basados en los periodos que deseamos explorar, por ejemplo hora del día,o día de la semana, o mes del año.

library(lubridate)
UKgrid_df$hour <- hour(UKgrid_df$time)
UKgrid_df$weekday <- wday(UKgrid_df$time, label = TRUE, abbr = TRUE)
UKgrid_df$month <- factor(month.abb[month(UKgrid_df$time)], levels =   month.abb)
head(UKgrid_df)
##                  time UKgrid hour weekday month
## 1 2005-04-01 00:00:00  65080    0     Fri   Apr
## 2 2005-04-01 01:00:00  68207    1     Fri   Apr
## 3 2005-04-01 02:00:00  69172    2     Fri   Apr
## 4 2005-04-01 03:00:00  66769    3     Fri   Apr
## 5 2005-04-01 04:00:00  64660    4     Fri   Apr
## 6 2005-04-01 05:00:00  65209    5     Fri   Apr

Vamos a empezar las exploraciones analizando el ciclo horario.

 UKgrid_hourly <- UKgrid_df %>%
    dplyr::group_by(hour) %>%
    dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE), sd = sd(UKgrid, na.rm
= TRUE))
str(UKgrid_hourly)
## tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
##  $ hour: int [1:24] 0 1 2 3 4 5 6 7 8 9 ...
##  $ mean: num [1:24] 55622 54343 52920 51546 50509 ...
##  $ sd  : num [1:24] 9259 9482 9485 9165 8557 ...
UKgrid_hourly
## # A tibble: 24 × 3
##     hour   mean     sd
##    <int>  <dbl>  <dbl>
##  1     0 55622.  9259.
##  2     1 54343.  9482.
##  3     2 52920.  9485.
##  4     3 51546.  9165.
##  5     4 50509.  8557.
##  6     5 51747.  8657.
##  7     6 59048. 10489.
##  8     7 68627. 13347.
##  9     8 73581. 13003.
## 10     9 76320. 12557.
## # … with 14 more rows
## # ℹ Use `print(n = ...)` to see more rows

Vamos ahora a hacer la gráfica de la media y la desviación estándar con base en las horas. Note que esas gráficas en escalas diferentes, así que hay que usar un gráfico especial.

require(plotly)
 plot_ly(UKgrid_hourly) %>%
      add_lines(x = ~ hour, y = ~ mean, name = "Media") %>%
      add_lines(x = ~ hour, y = ~ sd, name = "Desviación Estándar", yaxis =
   "y2",
                line = list(color = "red", dash = "dash", width = 3)) %>%
      layout(
        title = "La demanda nacional de electricidad - Promedio horario vs. Desviación Estándar",
        yaxis = list(title = "Media"),
        yaxis2 = list(overlaying = "y",
                      side = "right",
                      title = "Desviación Estándar"
        ),
        xaxis = list(title="Hora del Día"),
        legend = list(x = 0.05, y = 0.9),
        margin = list(l = 50, r = 50)
)

Qué podemos destacar del gráfico y de las estadísticas descriptivas?

  • Hay poca demanda durante la noche (entre la medianoche y las 6 a.m.) y una alta demanda entre las horas de la mañana y la tarde.

  • Existe una fuerte correlación entre la demanda promedio y su desviación estándar.

  • La relativamente baja desviación estándar de la demanda promedio durante la noche podría indicar que existe un fuerte efecto subestacional durante esas horas además de la estacionalidad horaria. Esto debería tener sentido, ya que son horas de sueño normales y, por lo tanto, en promedio, la demanda es razonablemente la misma durante los días de semana.

  • Por otro lado, la alta desviación estándar a lo largo de las horas de alta demanda podría indicar que la demanda se distribuye de manera diferente en diferentes vistas de periodicidad (como día de la semana o mes del año).

Vamos a explorar el último punto, viendo la demanda a la madrugada(3 a.m) y empezando el día(9 a.m) con respecto al día de la semana.

   UKgrid_weekday <- UKgrid_df %>%
      dplyr::filter(hour == 3 | hour == 9) %>%
    dplyr::group_by(hour, weekday) %>%
    dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
                     sd = sd(UKgrid, na.rm = TRUE))
## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.
UKgrid_weekday$hour <- factor(UKgrid_weekday$hour)
    plot_ly(data = UKgrid_weekday, x = ~ weekday, y = ~ mean, type =
   "bar",color = ~ hour) %>%
      layout(title = "The Hourly Average Demand by Weekday",
             yaxis = list(title = "Mean", range = c(30000, 75000)),
             xaxis = list(title = "Weekday"))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

En el gráfico de barras anterior podemos ver que la demanda de electricidad a las 3 a.m. es relativamente estable durante todos los días de la semana, con una ligera diferencia entre el promedio durante los días laborables y los días del fin de semana (alrededor de un 2% diferente). Por otro lado, existe una diferencia entre la demanda del día laborable y el fin de semana a las 9 a.m. (es decir, la demanda del lunes es en promedio un 28% superior a la del domingo). Como era de esperar, esos resultados se alinearon con nuestras expectativas anteriores.

Ahora podemos aprovechar esos conocimientos para examinar si existe un patrón estacional mensual en la serie. Ahora seleccionaremos las mismas horas (3 a.m. y 9 a.m.); sin embargo, esta vez agruparemos estos datos por mes (en lugar de días laborables):

 UKgrid_month <- UKgrid_df %>%
      dplyr::filter(hour == 3 | hour == 9) %>%
    dplyr::group_by(hour, month) %>%
    dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
                     sd = sd(UKgrid, na.rm = TRUE))
## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.
UKgrid_month$hour <- factor(UKgrid_month$hour)
    plot_ly(data = UKgrid_month, x = ~ month, y = ~ mean, type = "bar",color =
   ~ hour) %>%
      layout(title = "The Hourly Average Demand by Weekday",
             yaxis = list(title = "Mean", range = c(30000, 75000)),
             xaxis = list(title = "Month"))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Podemos ver en el gráfico de barras del resumen de agregación mensual que, en promedio, la demanda durante la noche (3 a.m.) y la mañana (9 a.m.) varía a lo largo de los meses del año. Además, hay un cambio significativo en la demanda durante la noche en comparación con la agregación entre semana. La variación de la serie de mes a mes indica la existencia de estacionalidad mensual en la serie.

Usando gráficos de densidades y Boxplot para explorar la estacionalidad.

Otro enfoque para analizar patrones estacionales en datos de series de tiempo es trazar la distribución de las unidades de frecuencia mediante el uso de histogramas o gráficos de densidad. Esto nos permitirá examinar si cada unidad de frecuencia tiene una distribución única que puede distinguirla del resto de unidades. Para esto usaremos el paquete ggplot2. Vamos empezar con la serie USgas.

require(timetk)
UKgrid_tbl
## # A tibble: 127,296 × 2
##    TIMESTAMP              ND
##    <dttm>              <int>
##  1 2005-04-01 00:00:00 65080
##  2 2005-04-01 01:00:00 68207
##  3 2005-04-01 02:00:00 69172
##  4 2005-04-01 03:00:00 66769
##  5 2005-04-01 04:00:00 64660
##  6 2005-04-01 05:00:00 65209
##  7 2005-04-01 06:00:00 72376
##  8 2005-04-01 07:00:00 82679
##  9 2005-04-01 08:00:00 88706
## 10 2005-04-01 09:00:00 91108
## # … with 127,286 more rows
## # ℹ Use `print(n = ...)` to see more rows
UKgrid_tbl%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = ND,.feature_set = c("hour","wday.lbl","month.lbl"))
## Warning: Removed 42 rows containing non-finite values (stat_boxplot).
UKgrid_tstbl%>%gg_subseries(ND,period=24)

 library(ggplot2)
    ggplot(USgas_df, aes(x = USgas)) +
      geom_density(aes(fill = month)) +
      ggtitle("USgas - Kernel Density Estimates by Month") +
      facet_grid(rows = vars(as.factor(month)))

Podemos ver algunos indicios de un patrón estacional en la serie, ya que las gráficas de densidad no se superponen entre sí (con la excepción de algunos meses consecutivos, como mayo y junio). Además, podemos ver que, durante algunos meses, la forma de las distribuciones es más plana con colas largas (principalmente durante los meses de invierno, noviembre, diciembre y enero). Sin embargo, no olvidemos el efecto de la tendencia o el crecimiento de un año a otro (como sabemos del capítulo anterior, la serie de gas de Estados Unidos tuvo una tendencia lineal desde el año 2010) ya que no la eliminamos de la serie. Repitamos este proceso; esta vez quitando la tendencia de la serie USgas antes de graficarla. Vamos a eliminarle la tendencia, posteriormente los explicaremos.

USgas_df$USgas_detrend <- USgas_df$USgas - decompose(USgas)$trend
    ggplot(USgas_df, aes(x = USgas_detrend)) +
      geom_density(aes(fill = month)) +
      ggtitle("USgas - Estimación de la densidad vía Kernel por mes") +
      facet_grid(rows = vars(as.factor(month)))
## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
## Warning: Removed 12 rows containing non-finite values (stat_density).

Podemos ver que hay un comportamiento similar pero ahora las colas de las densidades estimadas son mas cortas.Pueden hacer lo mismo pero en vez de quitar la tendencia, hacemos una diferenciación.

En el caso de que la distribución de la mayoría de las unidades de frecuencia sea plana con una cola larga, podría ser una indicación de múltiples patrones estacionales en la serie. Regresemos a la serie UKgrid y tracemos las gráficas de densidad de 24 horas:

UKgrid_df$hour <- as.factor(UKgrid_df$hour)
    ggplot(UKgrid_df, aes(x = UKgrid)) +
      geom_density(aes(fill = hour)) +
      ggtitle("UKgrid - Kernel Density Estimates by Hour of the day") +
      facet_grid(rows = vars(as.factor(hour)))
## Warning: Removed 14 rows containing non-finite values (stat_density).

Como observamos anteriormente con las tablas de resúmenes estadísticos, la distribución de la demanda neta de electricidad durante la noche es relativamente estable (de ahí la distribución no plana con colas cortas en contraposición a la distribución plana con cola larga durante el día). Si ahora hacemos un subconjunto con una de las horas durante el día y trazamos su distribución por el día de la semana, deberíamos esperar una superposición durante la noche y poder distinguir entre la distribución durante los días de la semana y el fin de semana, en contraposición a solo el día de la semana.

Por ejemplo, la siguiente gráfica representa la distribución de la demanda a las 9 a.m. a lo largo de los días de la semana. Puede ver que la distribución durante los días de la semana se distingue de la del fin de semana:

UKgrid_df$weekday <- as.factor(UKgrid_df$weekday)
    UKgrid_df %>% dplyr::filter(hour == 0) %>%
    ggplot(aes(x = UKgrid)) +
      geom_density(aes(fill = as.factor(weekday))) +
      ggtitle("UKgrid - Kernel Density Estimates by Hour of the day") +
      facet_grid(rows = vars(as.factor(weekday)))

Más herramientas del análisis de la estacionalidad

Vamos a usar el paquete forecast. Note que extraemos las subseries de los años. Lo cual demuestra un fuerte patrón estacional mensual.

library(forecast)
   ggseasonplot(USgas,year.labels=TRUE,continuous=TRUE)

ggseasonplot(USgas,  polar = TRUE)

Note que con el paquete TSstudio se puede hacer algo análogo.

ts_seasonal(USgas,type ="normal")
ts_seasonal(USgas, type = "cycle")
ts_seasonal(USgas, type = "box")
ts_seasonal(USgas, type = "all")

ts_seasonal con type=“cycle” añade un orden cronológico, en este caso por mes. Esto puede permitir la identificación de un patrón estacional sin tener que quitar la tendencia.

Cuando en el argumento type=“box”, elabora un gráfico de caja por unidad de frecuencia.

Finalmente cuando en el argumento type=“all”, elabora todos los gráficos anteriores.

El mapa de calor ya lo vimos anteriormente.

ts_heatmap(USgas, color = "Reds")

Gráficos basados en cuantiles también son de utilidad.De forma predeterminada, la función devuelve un gráfico de cuantiles de las unidades de frecuencia de la serie, donde la línea media representa la mediana y las líneas inferior y superior representan los percentiles 25 y 75.

ts_quantile(UKgrid)
## Warning in ts_quantile(UKgrid): The input object is a multiple time series
## object, by defualt will use only the first series as an input

El argumento period permite examinar si los patrones estacionales de la serie están cambiando cuando se usa un subconjunto de tiempo diferente. Esto le permite examinar si la serie tiene patrones estacionales adicionales. Por ejemplo, podemos trazar el ciclo de 24 horas de la serie UKgrid por día de la semana estableciendo el argumento de período en días de la semana:

ts_quantile(UKgrid, period = "weekdays", n = 2)
## Warning in ts_quantile(UKgrid, period = "weekdays", n = 2): The input object is
## a multiple time series object, by defualt will use only the first series as an
## input

Como vimos anteriormente con las gráficas de densidad, la demanda de electricidad durante el día es relativamente mayor durante los días de semana en comparación con los fines de semana. De la misma manera, puede trazar el ciclo de 24 horas por mes:

 ts_quantile(UKgrid, period = "monthly", n = 2)
## Warning in ts_quantile(UKgrid, period = "monthly", n = 2): The input object is
## a multiple time series object, by defualt will use only the first series as an
## input

La principal ventaja de las gráficas de cuantiles de múltiples períodos (por ejemplo, días de la semana, meses, etc.) sobre la gráfica de densidad que usamos anteriormente es que la primera representa todas las unidades de frecuencia (y, en el caso de UKgrid serie, con una frecuencia de 24 horas), mientras que el segundo representa una sola unidad de frecuencia (por ejemplo, la densidad de las 9 am durante los días de semana).

Usando Regresión para Descubrir un cíclo

Vamos a considerar la serie \[x_t =A\cos(2\pi\omega t+\varphi)+w_t,\] y una simulación de tamaño \(T=500\) para generar dos series, en donde las varianzas de los ruidos son diferentes en cada caso. Más específicamente \(\omega=\frac{1}{50}, A=2, \varphi=0.6\pi\) y \(\sigma_w=1,5.\)

cs = 2*cos(2*pi*1:500/50 + .6*pi);  w = rnorm(500,0,1)
par(mfrow=c(3,1), mar=c(3,2,2,1), cex.main=1.5)
plot.ts(cs, main=expression(2*cos(2*pi*t/50+.6*pi)))
plot.ts(cs+w, main=expression(2*cos(2*pi*t/50+.6*pi) + N(0,1)))
plot.ts(cs+5*w, main=expression(2*cos(2*pi*t/50+.6*pi) + N(0,25)))

Note que cuando la varianza del ruido es mas grande, es mas complicado ver el cíclo.

Vamos a re-escribir el proceso anterior usando la relación trigonométrica \(\cos(\alpha \pm \beta) = \cos(\alpha) \cos(\beta) \mp \sin(\alpha) \sin(\beta),\) con lo cual:

\[A \cos(2\pi\omega t + \varphi) = \beta_1 \cos(2\pi\omega t) + \beta_2 \sin(2\pi\omega t)\] donde \(\beta_1=A\cos(\varphi)\) y \(\beta_2=-A\sin(\varphi)\), así que el proceso queda escrito de la siguiente manera:

\[x_t = \beta_1 \cos(2\pi t/50) + \beta_2 \sin(2\pi t/50) + w_t.\] asimiendo que \(\omega=1/50\) es conocido. La transformación nos permitió convertir la regresión no lineal(en los parámetros) en una lineal. Así podemos estimar \(\beta_1\) y \(\beta_2\) usando una regresión lineal. Note también que si tomamos \(\beta_1^2+\beta_2^2\), entonces tenemos \(A^2=\beta_1^2+\beta_2^2\), así que, un estimador de la amplitud al cuadrado es \(\hat{A}^2=\hat{\beta}_1^2+\hat{\beta}_2^2.\) Note que los valores verdaderos de \(\beta_1=\) y \(\beta_2=\)

beta1teor=2*cos(.6*pi)
beta2teor=-2*sin(.6*pi)
set.seed(90210) # fijamos una semilla para que nos dé el mismo resultado
x = 2*cos(2*pi*1:500/10 + .6*pi) + rnorm(500,0,5)
z1 = cos(2*pi*1:500/50)
z2 = sin(2*pi*1:500/50)
summary(fit <- lm(x~0+z1+z2)) # sin intercepto.
## 
## Call:
## lm(formula = x ~ 0 + z1 + z2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5372  -4.0634  -0.2505   3.3048  16.1621 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)
## z1 -0.12622    0.34339  -0.368    0.713
## z2 -0.09283    0.34339  -0.270    0.787
## 
## Residual standard error: 5.429 on 498 degrees of freedom
## Multiple R-squared:  0.0004179,  Adjusted R-squared:  -0.003597 
## F-statistic: 0.1041 on 2 and 498 DF,  p-value: 0.9012
beta1teor
## [1] -0.618034
beta2teor
## [1] -1.902113
par(mar = c(2,2,2,2))
par(mfrow=c(2,1))
plot.ts(x)
plot.ts(x, col=8, ylab=expression(hat(x)))
  lines(fitted(fit), col=2)

El periodograma

El periodograma es una herramienta que permite detectar el valor de la frecuencia \(\omega\) que en general es desconocida.

Exploración de Múltiples Cíclos y el Periodograma

Vamos a explorar una herremienta que sirve para describir funciones periódicas generales. Fourier mostró que toda función periódica puede representarse como una suma de funciones sinusoidales de distinta amplitud y frecuencia. La idea entonces es generalizar el análisis anterior para un ciclo a suma de funciones armónicas con distintas frecuencias.

Dada la serie de longitud \(T\), se denomina periodo básico o de Fourier, a las fracciones exactas del tamaño muestral, es decir \[s_j=\frac{T}{j},\ \ j=1,2,\cdots,T/2.\]

El valor máximo de periodo básico se obtiene para \(j=1\) y es \(T\), es decir, observamos la onda una sola vez. El mínimo se obtiene para \(j=T/2\) y es 2 ya que no se puede observar periodos que duren menos de 2 observaciones. Como antes, suele usarse las frecuencias en lugar de los periodos para el ajuste de los ciclos, así que las frecuencias básicas o de Fourier se define como \[f_j=\frac{j}{T}\ \ j=1,2,\cdots,T/2.\]

Así, \(1/T\leq f_j\leq 1/2\), y el valor máximo sería \(f_j=0.5\) y se conoce como la frecuencia Nyquist. Con esto, una serie de tiempo \(Z_t\) que es periódica se puede representar mediante

\[Z_t=\mu+\sum_{j=1}^{T/2}A_j\sin(\omega_j t)+ \sum_{j=1}^{T/2}B_j\cos(\omega_j t)+a_t.\]

Note que éste modelo tiene tantos parámetros como observaciones, con lo cual hay que buscar un procedimiento para seleccionar las frecuencias que deben ser incluidas para explicar la evolución de la serie. Para esto se utilizará la herramienta conocida como el periodograma. Se puede verificar que la contribución de una onda a la varianza es dada por \(\frac{\hat{R}^2}{2}\), así que ondas asociadas con amplitudes grandes serán importantes en la explicación de la variabilidad de la serie, mientras que ondas asociadas con amplitudes bajas serán poco importantes. De forma análoga a como se mostró en la sección anterior, los estimadores de los coeficientes \(A_j\) y \(B_j\) están dados por \[\hat{A}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\sen(\omega_j t)\]
\[\hat{B}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\cos(\omega_j t)\] y \[\hat{R}_j=\hat{A}_j^2+\hat{B}_j^2,\] donde \(\omega_j=2\pi f_j\) y \(f=\frac{j}{T}\), con lo cual \[TS_z^2=\frac{T}{2}\sum_{j=1}^{T/2}\hat{R}_j^2.\]

Se le da el nombre de periodograma a la representación de cada \(\frac{T\hat{R}_j^2}{2}\) en función de la frecuencia \(\omega_j\) o \(f_j\). Así \[\begin{equation} \label{periodograma} I(f_j)=\frac{T\hat{R}_j^2}{2},\ \ \text{con}\ 1/T\leq f_j\leq 1/2 \end{equation}\] y \[\bar{I}=\sum_{j=1}^{T/2}\hat{R}_j^2=S_z^2.\]

Nos enfocaremos en las frecuencias básicas únicamente, lo cual no es restrictivo si \(T\) es muy grande, puesto que el número de frecuencias básicas será muy grande y siempre existirá alguna frecuencia básica muy próxima a la que puede interesarnos. El periodograma puede verse como una herramienta para la detección de posibles ciclos(ocultos) deterministas en una serie temporal. Por ejemplo, en una serie mensual estacional de periodo \(s=12\), esperamos encontrar un valor alto de periodograma para \(f=1/12\), pero también para \(f=j/12\), es decir, \(1/6,1/4,1/3\) que son armónicos del periodo estacional. Por otro lado, la serie puede tener otros ciclos no necesariamente ligado al periodo estacional, siendo el periodograma una buena herramienta para detectar estos posibles componentes. \

Podemos considerar la amplitud calculada para la frecuencia \(f_j\) como un promedio de las amplitudes existentes en las frecuencias situadas en el intervalo \(f_j\pm \frac{1}{2}T\). Con esto, se puede obtener el periodograma suavizado construyendo rectángulos con centro \(f_j\), base igual a \(1/T\) y alturas \(I(f_j)=T\hat{R}_j/2\). Otra forma de suavizar el periodograma es \[I(f)=\sum_{f_i-q}^{fi+q}p_iI(f_i)\ \ 0\leq f \leq 0.5,\]

donde \(q\) representa la ventana usada y los \(p_i\) son los pesos simétricos y no modifica el valor de la varianza \(S^2_z\) que es área bajo la curva.

Vamos a retomar el ejemplo anterior donde simulamos una serie de tamaño \(T=500\) con frecuencia \(\omega=1/50\).

plot.ts(x)

spectrum(x,log='no')
abline(v=1/10, lty=2,col="red")

###Otra forma de computarlo es vía la transformada de Fourier

Px = Mod(fft(x))^2
Freq=0:499/500
plot(Freq, Px, type='l')

u = which.max(Px[1:250])
sprintf("El valor de la frecuencia donde se máximiza el periodograma es %s",Freq[u])
## [1] "El valor de la frecuencia donde se máximiza el periodograma es 0.1"

Note que en cualquiera de las dos formas de obtener el periodograma, hemos descubierto la frecuencia \(\omega=1/50\) como frecuencia príncipal.

En ocasiones es necesario suavizar el periodograma porque encontramos varios picos que en verdad no son valores grandes.

spectrum(x,log='no',span=5)

#span vector of odd integers giving the widths of modified Daniell smoothers to be used to smooth the periodogram.

spectrum(x,log='no',span=c(5,5))

spectrum(x,log='no',span=c(2,2))

Qué sucede si no consideramos un ciclo determinístico a través de funciones sinusoidales sino que lo consideramos a través de variables Dummy? Es decir por ejemplo una componente estacional de periodo 4, modelado a través de variables dummy. \[Y_t=\delta_1\gamma_{1,t}+\delta_2\gamma_{2,t}+\delta_3\gamma_{3,t}+\delta_4\gamma_{4,t}+\epsilon_{t}\] donde las variables dummy \(\gamma_{i,t}, i=1,2,3,4\) indican si se está en el trimestre \(i\). Puede en periodograma detactar este tipo de estacionalidades obtenidas a través de dummy estacionales?

library(dplyr)
library(rstanarm)
## Loading required package: Rcpp
## This is rstanarm version 2.21.1
## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
## - For execution on a local, multicore CPU with excess RAM we recommend calling
##   options(mc.cores = parallel::detectCores())
## 
## Attaching package: 'rstanarm'
## The following object is masked from 'package:car':
## 
##     logit
## The following objects are masked from 'package:VGAM':
## 
##     cauchy, dirichlet, exponential, laplace, logit
library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
library(ggplot2)
library(bayesplot)
## This is bayesplot version 1.9.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
s=4
deltas=seq(1:s)

####Ejemplo Fácil Dummy+ruido

T=200*s
epsilon=rnorm(T,0,1)
l=diag(s)
X=matrix(rep(t(l),T/s),ncol=ncol(l),byrow=TRUE)
Y=X%*%deltas+epsilon
plot(as.ts(Y))

acf(as.ts(Y))

#simul=data.frame(resp=Y,diseno=X)
#salida_dummy_ruido=stan_glm(resp~. -1,data=simul)
PeriodgramaY=spectrum(Y,log='no')

ubicacionY=which.max(PeriodgramaY$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para la serie es: %s",PeriodgramaY$freq[ubicacionY])
## [1] "El valor de la frecuencia donde se máximiza el periodograma para la serie es: 0.25"
sprintf("El periodo correspondiente es aporximadamente: %s",1/PeriodgramaY$freq[ubicacionY])
## [1] "El periodo correspondiente es aporximadamente: 4"

Veamos ahora con una serie de tiempo real. Usaremos la serie SOi y Recruitment para ilustrar su uso. Recordemos que ambas series son mensuales, así que las frecuencias que serán dibujadas en el eje x vendrán en múltiplos de \(\frac{1}{12}\). Usaremos ahora la función mvspec del paquete astsa.

library(astsa)
data(soi)
soi.per = astsa::mvspec(soi, log="no")
ubicacionsoi=which.max(soi.per$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para SOI es: %s",soi.per$freq[ubicacionsoi])
## [1] "El valor de la frecuencia donde se máximiza el periodograma para SOI es: 1"
n_soi <- length(soi.per$spec)
valor_seg_per=sort(soi.per$spec,partial=n_soi-1)[n_soi-1]
ubica_segundo_soi=which(soi.per$spec==valor_seg_per)

sprintf("El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para SOI es: %s",soi.per$freq[ubica_segundo_soi])
## [1] "El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para SOI es: 0.2"
abline(v=soi.per$freq[ubicacionsoi], lty=2,col="blue")
abline(v=soi.per$freq[ubica_segundo_soi], lty=2,col="blue")

Note que la primera frecuencia mas alta se encuentra en \(\omega=1\cdot \frac{1}{12}=\frac{1}{12}\), el cual corresponde al obvio ciclo anual (periodo de 12 meses). La segunda corresponde a \(\omega=\frac{1}{5} \cdot \frac{1}{12}=\frac{1}{60}\), el cual corresponde a un posible ciclo de periodo \(1/(1/60)= 60\) meses, es decir 5 años, lo cual puede ser debido al fenómeno de El Niño.

data("rec")
rec.per = astsa::mvspec(rec, log="no")
ubicacionrec=which.max(rec.per$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para REC es: %s",rec.per$freq[ubicacionrec])
## [1] "El valor de la frecuencia donde se máximiza el periodograma para REC es: 1"
n_rec <- length(rec.per$spec)
valor_seg_per_rec=sort(rec.per$spec,partial=n_rec-1)[n_rec-1]
ubica_segundo_rec=which(rec.per$spec==valor_seg_per_rec)

sprintf("El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para REC es: %s",rec.per$freq[ubica_segundo_rec])
## [1] "El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para REC es: 0.25"
abline(v=rec.per$freq[ubicacionrec], lty=2,col="blue")
abline(v=rec.per$freq[ubica_segundo_rec], lty=2,col="blue")

Note que la primera frecuencia mas alta se encuentra en \(\omega=1\cdot \frac{1}{12}=\frac{1}{12}\), el cual corresponde al obvio ciclo anual (periodo de 12 meses), igual que para el casi del indice soi. La segunda corresponde a \(\omega=\frac{1}{4} \cdot \frac{1}{12}=\frac{1}{48}\), el cual corresponde a un posible ciclo de periodo \(1/(1/48)= 48\) meses, es decir 4 años. lo cual puede ser debido al fenómeno de El Niño. Esta actividad de banda ancha sugiere que el posible ciclo de El Niño es irregular, pero tiende a rondar los cuatro años en promedio.

Tarea : Hacer lo mismo pero ahora usar la serie de pasajeros diferenciada, al igual que la serie UKgrid, mas específicamente la serie UKgrid_xts y la serie taylor del paquete forecast.

library(sarima) 
## 
## Attaching package: 'sarima'
## The following object is masked from 'package:rstanarm':
## 
##     se
## The following object is masked from 'package:astsa':
## 
##     sarima
## The following object is masked from 'package:stats':
## 
##     spectrum
x <- ts(sarima::sim_sarima(n=144, model = list(siorder=1,
                                      nseasons=12, sigma2 = 1)),frequency = 12 )#(1-B)^{12} X_t = e_t
plot(x)

acf(x)

monthplot(x)

ts_seasonal(x, type = "all")
spectrum(x,log='no',span=5)
## Estimated spectral density
##     series: x
##     method: Smoothed Periodogram
##     nseasons: 12
##     frequency range: (-6,6]
## 
##     sort method for the table: decreasing magnitudes
## 
##               freq      spec    % Total     Cum. %
##     [1,] 1.0000000 0.5488446 0.09121551 0.09121551
##     [2,] 0.9166667 0.5259173 0.08740510 0.17862061
##     [3,] 1.0833333 0.5008903 0.08324572 0.26186633
##     [4,] 0.8333333 0.3152458 0.05239244 0.31425877
##     [5,] 2.0833333 0.3137395 0.05214210 0.36640087
##     [6,] 2.0000000 0.3113759 0.05174928 0.41815014

astsa::mvspec(x, log="no")

library(sarima) 
x <- ts(sim_sarima(n=288,model=list(sar=0.3, nseasons=12, sigma2 = 1)), frequency = 12 )#SAR(1)
plot(x)

acf(x)

monthplot(x)

ts_seasonal(x, type = "all")
spectrum(x,log='no',span=5)
## Estimated spectral density
##     series: x
##     method: Smoothed Periodogram
##     nseasons: 12
##     frequency range: (-6,6]
## 
##     sort method for the table: decreasing magnitudes
## 
##              freq      spec    % Total     Cum. %
##     [1,] 2.000000 0.2843993 0.02370651 0.02370651
##     [2,] 2.041667 0.2799563 0.02333616 0.04704267
##     [3,] 4.166667 0.2643718 0.02203709 0.06907975
##     [4,] 4.083333 0.2621709 0.02185363 0.09093338
##     [5,] 4.125000 0.2549984 0.02125575 0.11218914
##     [6,] 1.958333 0.2365108 0.01971469 0.13190383

astsa::mvspec(x, log="no")

raices=polyroot(c(1,rep(0,11),-0.3))
raices
##  [1]  0.5527684+0.9574230i -0.9574230+0.5527684i -0.5527684-0.9574230i
##  [4]  0.9574230-0.5527684i  0.0000000+1.1055369i -1.1055369-0.0000000i
##  [7]  0.0000000-1.1055369i  1.1055369+0.0000000i -0.5527684+0.9574230i
## [10] -0.9574230-0.5527684i  0.5527684-0.9574230i  0.9574230+0.5527684i
conj_raices=Conj(raices)
raices[12]*conj_raices[12]
## [1] 1.222212+0i

Suavizamiento en el contexto de Series de Tiempo

Promedio Móvil

El promedio móvil es un método es útil para descubrir ciertos rasgos en una serie de tiempo, como tendencias a largo plazo y componentes estacionales. En particular, si \(x_t\) representa las observaciones, entonces una forma de predecir predecir o estimar la tendencia de la serie es:

\[m_t=\sum_{j=-k}^{k}a_jx_{t-j},\] donde si\(a_j=a_{-j}\ge0\) y \(\sum_{j=-k}^{k}a_j=1\) se conoce como el promedio móvil simétrico de los datos.

# Filtro Boxcar
wgts = c(.5, rep(1,11), .5)/12   ###Los pesos de los filtros
soif = stats::filter(soi, sides=2, filter=wgts)
plot(soi)
lines(soif, lwd=2, col=4)

###Forma del Filtro
nwgts = c(rep(0,20), wgts, rep(0,20))
plot(nwgts, type="l", ylim = c(-.02,.1), xaxt='n', yaxt='n', ann=FALSE)

## Suavizamiento Kernel El suavizamiento kernel es un suavizador de promedio móvil que utiliza una función de ponderación, o kernel, para promediar las observaciones.Veamos ahora como queda el promedio móvil:

\[m_t=\sum_{i=1}^{n}w_i(t)x_i\] donde \[w_i(t)=K(\frac{t-i}{b})/\sum_{j=1}^{n}K(\frac{t-j}{b})\] son los pesos y \(K(\cdot)\) es una función kernel. Este estimador es llamado el Estimador de Nadaraya-Watson. Usaremos la función ksmooth

plot(soi)
lines(ksmooth(time(soi), soi, "normal", bandwidth=1), lwd=2, col=4)

#par(fig = c(.65, 1, .65, 1), new = TRUE) # the insert
gauss = function(x) { 1/sqrt(2*pi) * exp(-(x^2)/2) }
x = seq(from = -3, to = 3, by = 0.001)
plot(x, gauss(x), type ="l", ylim=c(-.02,.45), xaxt='n', yaxt='n', ann=FALSE)

## Lowess Otro enfoque para suavizar un gráfico de tiempo es la regresión del vecino más cercano. La técnica se basa en la regresión de k vecinos más cercanos, en la que uno usa solo los datos \(\{x_{t−k/2}, ..., x_t, ..., x_{t+k/2}\}\) para predecir \(x_t\) mediante regresión, y luego establece \(m_t = \hat{x}_t\). Primero, una cierta proporción de vecinos más cercanos a \(x_t\) se incluyen en un esquema de ponderación; los valores más cercanos a \(x_t\) en el tiempo obtienen más peso. Luego, se utiliza una regresión ponderada robusta para predecir \(x_t\) y obtener los valores suavizados \(mt\). Cuanto mayor sea la fracción de vecinos más cercanos incluidos, más suave será el ajuste. El R se usa la función lowess.

plot(soi)
lines(lowess(soi, f=.05), lwd=2, col=4)  # Ciclo El Niño 
lines(lowess(soi), lty=2, lwd=2, col=2)  # tendencia (con span por defecto)

Suavizamiento Splines

Una forma obvia de suavizar los datos sería ajustar una regresión polinomial en términos del tiempo. Por ejemplo, un polinomio cúbico tendría \(x_t = m_t + w_t\) donde \(m_t =\beta_0 + \beta_1t + \beta_2t^2 + \beta_3t^3\). Entonces podríamos ajustar \(m_t\) mediante mínimos cuadrados ordinarios.

Una extensión de la regresión polinomial es dividir primero el tiempo \(t = 1,. . . , n\), en k intervalos, \([t_0 = 1, t_1]\), \([t_1 + 1, t_2],\cdots,\) \([t_{k − 1} + 1, t_k = n]\); los valores \(t_0\), \(t_1\), …, \(t_k\) se llaman nodos. Luego, en cada intervalo, se ajusta una regresión polinomial, normalmente de orden 3, y esto se llama splines cúbicos. Un método relacionado es suavizar splines, que minimiza el compromiso entre el ajuste y el grado de suavidad dado por

\[\sum_{t=1}^{n}[x_t-m_t]^2+\lambda\int(m_t'')^2 dt,\] donde \(m_t\) es un spline cúbico con nodos en cada tiempo t y el grado de suavidad es controlado por \(\lambda>0.\) El parámetro de suavizado en R es controlado por el argumento spar de la función smooth.spline del paquete stats.

plot(soi)
  lines(smooth.spline(time(soi), soi, spar=.5), lwd=2, col=4)
  lines(smooth.spline(time(soi), soi, spar= 1), lty=2, lwd=2, col=2)

La diferencia entre uan y otra gráfica es es grado de suavizamiento.

<> Tarea: Hacer lo mismo para la serie Recruitment.

#```{r} #library(“readr”) #library(zoo)

#write_csv(data.frame(soi=as.matrix(soi),date=as.Date(as.yearmon(time(soi)))), file = “soi.csv”)

#write_csv(data.frame(rec=as.matrix(rec),date=as.Date(as.yearmon(time(rec)))), file = “rec.csv”) #```